RPYFMM: Parallel Adaptive Fast Multipole Method for Rotne-Prager-Yamakawa Tensor in Biomolecular Hydrodynamics Simulations.

RPYFMM is a software package for the efficient evaluation of the potential field governed by the Rotne-Prager-Yamakawa (RPY) tensor interactions in biomolecular hydrodynamics simulations. In our algorithm, the RPY tensor is decomposed as a linear combination of four Laplace interactions, each of which is evaluated using the adaptive fast multipole method (FMM) [1] where the exponential expansions are applied to diagonalize the multipole-to-local translation operators. RPYFMM offers a unified execution on both shared and distributed memory computers by leveraging the DASHMM library [2, 3]. Preliminary numerical results show that the interactions for a molecular system of 15 million particles (beads) can be computed within one second on a Cray XC30 cluster using 12, 288 cores, while achieving approximately 54% strong-scaling efficiency.

Enzyme localization, crowding, and buffers collectively modulate diffusion-influenced signal transduction: Insights from continuum diffusion modeling.

Biochemical reaction networks consisting of coupled enzymes connect substrate signaling events with biological function. Substrates involved in these reactions can be strongly influenced by diffusion "barriers" arising from impenetrable cellular structures and macromolecules, as well as interactions with biomolecules, especially within crowded environments. For diffusion-influenced reactions, the spatial organization of diffusion barriers arising from intracellular structures, non-specific crowders, and specific-binders (buffers) strongly controls the temporal and spatial reaction kinetics. In this study, we use two prototypical biochemical reactions, a Goodwin oscillator, and a reaction with a periodic source/sink term to examine how a diffusion barrier that partitions substrates controls reaction behavior. Namely, we examine how conditions representative of a densely packed cytosol, including reduced accessible volume fraction, non-specific interactions, and buffers, impede diffusion over nanometer length-scales. We find that diffusion barriers can modulate the frequencies and amplitudes of coupled diffusion-influenced reaction networks, as well as give rise to "compartments" of decoupled reactant populations. These effects appear to be intensified in the presence of buffers localized to the diffusion barrier. These findings have strong implications for the role of the cellular environment in tuning the dynamics of signaling pathways.

Predicting the influence of long-range molecular interactions on macroscopic-scale diffusion by homogenization of the Smoluchowski equation.

The macroscopic diffusion constant for a charged diffuser is in part dependent on (1) the volume excluded by solute "obstacles" and (2) long-range interactions between those obstacles and the diffuser. Increasing excluded volume reduces transport of the diffuser, while long-range interactions can either increase or decrease diffusivity, depending on the nature of the potential. We previously demonstrated [P. M. Kekenes-Huskey et al., Biophys. J. 105, 2130 (2013)] using homogenization theory that the configuration of molecular-scale obstacles can both hinder diffusion and induce diffusional anisotropy for small ions. As the density of molecular obstacles increases, van der Waals (vdW) and electrostatic interactions between obstacle and a diffuser become significant and can strongly influence the latter's diffusivity, which was neglected in our original model. Here, we extend this methodology to include a fixed (time-independent) potential of mean force, through homogenization of the Smoluchowski equation. We consider the diffusion of ions in crowded, hydrophilic environments at physiological ionic strengths and find that electrostatic and vdW interactions can enhance or depress effective diffusion rates for attractive or repulsive forces, respectively. Additionally, we show that the observed diffusion rate may be reduced independent of non-specific electrostatic and vdW interactions by treating obstacles that exhibit specific binding interactions as "buffers" that absorb free diffusers. Finally, we demonstrate that effective diffusion rates are sensitive to distribution of surface charge on a globular protein, Troponin C, suggesting that the use of molecular structures with atomistic-scale resolution can account for electrostatic influences on substrate transport. This approach offers new insight into the influence of molecular-scale, long-range interactions on transport of charged species, particularly for diffusion-influenced signaling events occurring in crowded cellular environments.

Finite Element Estimation of Protein-Ligand Association Rates with Post-Encounter Effects: Applications to Calcium binding in Troponin C and SERCA.

We introduce a computational pipeline and suite of software tools for the approximation of diffusion-limited binding based on a recently developed theoretical framework. Our approach handles molecular geometries generated from high-resolution structural data and can account for active sites buried within the protein or behind gating mechanisms. Using tools from the FEniCS library and the APBS solver, we implement a numerical code for our method and study two Ca(2+)-binding proteins: Troponin C and the Sarcoplasmic Reticulum Ca(2+) ATPase (SERCA). We find that a combination of diffusional encounter and internal 'buried channel' descriptions provide superior descriptions of association rates, improving estimates by orders of magnitude.

Dewetting-controlled binding of ligands to hydrophobic pockets.

We report on a combined atomistic molecular dynamics simulation and implicit solvent analysis of a generic hydrophobic pocket-ligand (host-guest) system. The approaching ligand induces complex wetting-dewetting transitions in the weakly solvated pocket. The transitions lead to bimodal solvent fluctuations which govern magnitude and range of the pocket-ligand attraction. A recently developed implicit water model, based on the minimization of a geometric functional, captures the sensitive aqueous interface response to the concave-convex pocket-ligand configuration semiquantitatively.

A virtual screening study of the acetylcholine binding protein using a relaxed-complex approach.

The nicotinic acetylcholine receptor (nAChR) is a member of the ligand-gated ion channel family and is implicated in many neurological events. Yet, the receptor is difficult to target without high-resolution structures. In contrast, the structure of the acetylcholine binding protein (AChBP) has been solved to high resolution, and it serves as a surrogate structure of the extra-cellular domain in nAChR. Here we conduct a virtual screening study of the AChBP using the relaxed-complex method, which involves a combination of molecular dynamics simulations (to achieve receptor structures) and ligand docking. The library screened through comes from the National Cancer Institute, and its ligands show great potential for binding AChBP in various manners. These ligands mimic the known binders of AChBP; a significant subset docks well against all species of the protein and some distinguish between the various structures. These novel ligands could serve as potential pharmaceuticals in the AChBP/nAChR systems.

Finite element analysis of drug electrostatic diffusion: inhibition rate studies in N1 neuraminidase.

This article describes a numerical solution of the steady-state Poisson-Boltzmann-Smoluchowski (PBS) and Poisson-Nernst-Planck (PNP) equations to study diffusion in biomolecular systems. Specifically, finite element methods have been developed to calculate electrostatic interactions and ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to the wild-type and several mutated avian influenza neurominidase crystal structures. The calculated rates show very good agreement with recent experimental studies. Furthermore, these finite element methods require significantly fewer computational resources than existing particle-based Brownian dynamics methods and are robust for complicated geometries. The key finding of biological importance is that the electrostatic steering plays the important role in the drug binding process of the neurominidase.

Acetylcholinesterase: mechanisms of covalent inhibition of H447I mutant determined by computational analyses.

The reaction mechanisms of two inhibitor TFK(+) and TFK(0) binding to H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. TFK(+) binding to the H447I mutant may proceed with a different reaction mechanism from the wild-type. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK(0), however, multiple MD simulations on the TFK(0)/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Taken together our computational studies confirm that TFK(0) is almost inactive in the H447I mutant, and also provide detailed mechanistic insights into the experimental observations.

Optimizing the Poisson Dielectric Boundary with Explicit Solvent Forces and Energies: Lessons Learned with Atom-Centered Dielectric Functions.

Accurate implicit solvent models require parameters that have been optimized in a system- or atom-specific manner on the basis of experimental data or more rigorous explicit solvent simulations. Models based on the Poisson or Poisson-Boltzmann equation are particularly sensitive to the nature and location of the boundary which separates the low dielectric solute from the high dielectric solvent. Here, we present a novel method for optimizing the solute radii, which define the dielectric boundary, on the basis of forces and energies from explicit solvent simulations. We use this method to optimize radii for protein systems defined by AMBER ff99 partial charges and a spline-smoothed solute surface. The spline-smoothed surface is an atom-centered dielectric function that enables stable and efficient force calculations. We explore the relative performance of radii optimized with forces alone and those optimized with forces and energies. We show that our radii reproduce the explicit solvent forces and energies more accurately than four other parameter sets commonly used in conjunction with the AMBER force field, each of which has been appropriately scaled for spline-smoothed surfaces. Finally, we demonstrate that spline-smoothed surfaces show surprising accuracy for small, compact systems but may have limitations for highly solvated protein systems. The optimization method presented here is efficient and applicable to any system with explicit solvent parameters. It can be used to determine the optimal continuum parameters when experimental solvation energies are unavailable and the computational costs of explicit solvent charging free energies are prohibitive.

A minimal model for stabilization of biomolecules by hydrocarbon cross-linking.

Programmed cell death regulating protein motifs play an essential role in the development of an organism, its immune response, and disease-related cellular mechanisms. Among those motifs the BH3 domain of the BCL-2 family is found to be of crucial importance. Recent experiments showed how the isolated, otherwise unstructured BH3 peptide can be modified by a hydrocarbon linkage to regain function. We parametrized a reduced, dynamic model for the stability effects of such covalent cross-linking and confirmed that the model reproduces the reinforcement of the structural stability of the BH3 motif by cross-linking. We show that an analytically solvable model for thermostability around the native state is not capable of reproducing the stabilization effect. This points to the crucial importance of the peptide dynamics and the fluctuations neglected in the analytic model for the cross-linking system to function properly. This conclusion is supported by a thorough analysis of a simulated Go model. The resulting model is suitable for rational design of generic cross-linking systems in silicio.

Coupling hydrophobicity, dispersion, and electrostatics in continuum solvent models.

An implicit solvent model is presented that couples hydrophobic, dispersion, and electrostatic solvation energies by minimizing the system Gibbs free energy with respect to the solvent volume exclusion function. The solvent accessible surface is the output of the theory. The method is illustrated with the solvation of simple solutes on different length scales and captures the sensitivity of hydration to the particular form of the solute-solvent interactions in agreement with recent computer simulations.

Coupling nonpolar and polar solvation free energies in implicit solvent models.

Recent studies on the solvation of atomistic and nanoscale solutes indicate that a strong coupling exists between the hydrophobic, dispersion, and electrostatic contributions to the solvation free energy, a facet not considered in current implicit solvent models. We suggest a theoretical formalism which accounts for coupling by minimizing the Gibbs free energy of the solvent with respect to a solvent volume exclusion function. The resulting differential equation is similar to the Laplace-Young equation for the geometrical description of capillary interfaces but is extended to microscopic scales by explicitly considering curvature corrections as well as dispersion and electrostatic contributions. Unlike existing implicit solvent approaches, the solvent accessible surface is an output of our model. The presented formalism is illustrated on spherically or cylindrically symmetrical systems of neutral or charged solutes on different length scales. The results are in agreement with computer simulations and, most importantly, demonstrate that our method captures the strong sensitivity of solvent expulsion and dewetting to the particular form of the solvent-solute interactions.

E230Q mutation of the catalytic subunit of cAMP-dependent protein kinase affects local structure and the binding of peptide inhibitor.

The active site of the mammalian cAMP-dependent protein kinase catalytic subunit (C-subunit) has a cluster of nonconserved acidic residues-Glu127, Glu170, Glu203, Glu230, and Asp241-that are crucial for substrate recognition and binding. Studies have shown that the Glu230 to Gln mutant (E230Q) of the enzyme has physical properties similar to the wild-type enzyme and has decreased affinity for a short peptide substrate, Kemptide. However, recent experiments intended to crystallize ternary complex of the E230Q mutant with MgATP and protein kinase inhibitor (PKI) could only obtain crystals of the apo-enzyme of E230Q mutant. To deduce the possible mechanism that prevented ternary complex formation, we used the relaxed-complex method (Lin, J.-H., et al. J Am Chem Soc 2002, 24, 5632-5633) to study PKI binding to the E230Q mutant C-subunit. In the E230Q mutant, we observed local structural changes of the peptide binding site that correlated closely to the reduced PKI affinity. The structural changes occurred in the F-to-G helix loop and appeared to hinder PKI binding. Reduced electrostatic potential repulsion among Asp241 from the helix loop section and the other acidic residues in the peptide binding site appear to be responsible for the structural change.

Electrostatic properties of cowpea chlorotic mottle virus and cucumber mosaic virus capsids.

Electrostatic properties of cowpea chlorotic mottle virus (CCMV) and cucumber mosaic virus (CMV) were investigated using numerical solutions to the Poisson-Boltzmann equation. Experimentally, it has been shown that CCMV particles swell in the absence of divalent cations when the pH is raised from 5 to 7. CMV, although structurally homologous, does not undergo this transition. An analysis of the calculated electrostatic potential confirms that a strong electrostatic repulsion at the calcium-binding sites in the CCMV capsid is most likely the driving force for the capsid swelling process during the release of calcium. The binding interaction between the encapsulated genome material (RNA) inside of the capsid and the inner capsid shell is weakened during the swelling transition. This probably aids in the RNA release process, but it is unlikely that the RNA is released through capsid openings due to unfavorable electrostatic interaction between the RNA and capsid inner shell residues at these openings. Calculations of the calcium binding energies show that Ca(2+) can bind both to the native and swollen forms of the CCMV virion. Favorable binding to the swollen form suggests that Ca(2+) ions can induce the capsid contraction and stabilize the native form.

Computing the Amino Acid Specificity of Fluctuations in Biomolecular Systems.

We developed a new amino acid specific method for the computation of spatial fluctuations of proteins around their native structures. We show the consistency with experimental values and the increased performance in comparison to an established model, based on statistical estimates for a set of test proteins. We apply the new method to HIV-1 protease in its wild-type form and to a V82F-I84V mutant that shows resistance to protease inhibitors. We further show how the method can be successfully used to explain the molecular biophysics of drug resistance of the mutant.

Substrate concentration dependence of the diffusion-controlled steady-state rate constant.

The Smoluchowski approach to diffusion-controlled reactions is generalized to interacting substrate particles by including the osmotic pressure and hydrodynamic interactions of the nonideal particles in the Smoluchoswki equation within a local-density approximation. By solving the strictly linearized equation for the time-independent case with absorbing boundary conditions, we present an analytic expression for the diffusion-limited steady-state rate constant for small substrate concentrations in terms of an effective second virial coefficient B2*. Comparisons to Brownian dynamics simulations excluding hydrodynamic interactions show excellent agreement up to bulk number densities of B2*rho0 < approximately = 0.4 for hard sphere and repulsive Yukawa-like interactions between the substrates. Our study provides an alternative way to determine the second virial coefficient of interacting macromolecules experimentally by measuring their steady-state rate constant in diffusion-controlled reactions at low densities.

Native-state conformational dynamics of GART: a regulatory pH-dependent coil-helix transition examined by electrostatic calculations.

Glycinamide ribonucleotide transformylase (GART) undergoes a pH-dependent coil-helix transition with pK(a) approximately 7. An alpha-helix is formed at high pH spanning 8 residues of a 21-residue-long loop, comprising the segment Thr120-His121-Arg122-Gln123-Ala124-Leu125-Glu126-Asn127. To understand the electrostatic nature of this loop-helix, called the activation loop-helix, which leads to the formation and stability of the alpha-helix, pK(a) values of all ionizable residues of GART have been calculated, using Poisson-Boltzmann electrostatic calculations and crystallographic data. Crystallographic structures of high and low pH E70A GART have been used in our analysis. Low pK(a) values of 5.3, 5.3, 3.9, 1.7, and 4.7 have been calculated for five functionally important histidines, His108, His119, His121, His132, and His137, respectively, using the high pH E70A GART structure. Ten theoretical single and double mutants of the high pH E70A structure have been constructed to identify pairwise interactions of ionizable residues, which have aided in elucidating the multiplicity of electrostatic interactions of the activation loop-helix, and the impact of the activation helix on the catalytic site. Based on our pK(a) calculations and structural data, we propose that: (1) His121 forms a molecular switch for the coil-helix transition of the activation helix, depending on its protonation state; (2) a strong electrostatic interaction between His132 and His121 is observed, which can be of stabilizing or destabilizing nature for the activation helix, depending on the relative orientation and protonation states of the rings of His121 and His132; (3) electrostatic interactions involving His119 and Arg122 play a role in the stability of the activation helix; and (4) the activation helix contains the helix-promoting sequence Arg122-Gln123-Ala124-Leu125-Glu126, but its alignment relative to the N and C termini of the helix is not optimal, and is possibly of a destabilizing nature. Finally, we provide electrostatic evidence that the formation and closure of the activation helix create a hydrophobic environment for catalytic-site residue His108, to facilitate catalysis.

Proton transfer dynamics of GART: the pH-dependent catalytic mechanism examined by electrostatic calculations.

The enzyme glycinamide ribonucleotide transformylase (GART) catalyzes the transfer of a formyl group from formyl tetrahydrofolate (fTHF) to glycinamide ribonucleotide (GAR), a process that is pH-dependent with pK(a) of approximately 8. Experimental studies of pH-rate profiles of wild-type and site-directed mutants of GART have led to the proposal that His108, Asp144, and GAR are involved in catalysis, with His108 being an acid catalyst, while forming a salt bridge with Asp144, and GAR being a nucleophile to attack the formyl group of fTHF. This model implied a protonated histidine with pK(a) of 9.7 and a neutral GAR with pK(a) of 6.8. These proposed unusual pK(a)s have led us to investigate the electrostatic environment of the active site of GART. We have used Poisson-Boltzmann-based electrostatic methods to calculate the pK(a)s of all ionizable groups, using the crystallographic structure of a ternary complex of GART involving the pseudosubstrate 5-deaza-5,6,7,8-THF (5dTHF) and substrate GAR. Theoretical mutation and deletion analogs have been constructed to elucidate pairwise electrostatic interactions between key ionizable sites within the catalytic site. Also, a construct of a more realistic catalytic site including a reconstructed pseudocofactor with an attached formyl group, in an environment with optimal local van der Waals interactions (locally minimized) that imitates closely the catalytic reactants, has been used for pK(a) calculations. Strong electrostatic coupling among catalytic residues His108, Asp144, and substrate GAR was observed, which is extremely sensitive to the initial protonation and imidazole ring flip state of His108 and small structural changes. We show that a proton can be exchanged between GAR and His108, depending on their relative geometry and their distance to Asp144, and when the proton is attached on His108, catalysis could be possible. Using the formylated locally minimized construct of GART, a high pK(a) for His108 was calculated, indicating a protonated histidine, and a low pK(a) for GAR(NH(2)) was calculated, indicating that GAR is in neutral form. Our results are in qualitative agreement with the current mechanistic picture of the catalytic process of GART deduced from the experimental data, but they do not reproduce the absolute magnitude of the pK(a)s extracted from fits of k(cat)-pH profiles, possibly because the static time-averaged crystallographic structure does not describe adequately the dynamic nature of the catalytic site during binding and catalysis. In addition, a strong effect on the pK(a) of GAR(NH(2)) is produced by the theoretical mutations of His108Ala and Asp144Ala, which is not in agreement with the observed insensitivity of the pK(a) of GAR(NH(2)) modeled from the experimental data using similar mutations. Finally, we show that important three-way electrostatic interactions between highly conserved His137, with His108 and Asp144, are responsible for stabilizing the electrostatic microenvironment of the catalytic site. In conclusion, our data suggest that further detailed computational and experimental work is necessary.

Ordered water and ligand mobility in the HIV-1 integrase-5CITEP complex: a molecular dynamics study.

A 2 ns molecular dynamics simulation has been carried out for the HIV-1 integrase-5CITEP complex in order to understand the role of water in defining the ligand's binding mode and to address issues of binding site flexibility and ligand motion. Although the ligand retains considerable mobility within the active site, a structural water molecule bridging 5CITEP with Asp 64 and Asn 155 is identified in the simulation. Consideration of this water molecule could open a route to new HIV-1 integrase inhibitors.